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Abstract 

We point out that the local density approximation (LDA) of Oliva is 
an adaptation of the Thomas-Fermi method, and is a good approximation 
when e = hu/kT << 1. For the case of scattering length a > 0, the LDA 
leads to a quantitative result (14') easily checked by experiments. Critical 
remarks are made about the physics of the many body problem in terms 
of the scattering length a. 



Bose-Einstein condensation (BEC) of free particles was a great contribution 
to physics which provides one more example of the awesome and daring insight of 
Einstein. For some sixty years, experimental observation of BEC was considered 
hopeless. Now through a series of ingenious developments BEC has finally been 
observed [1-3] for trapped bosons. We would like to discuss here BEC of such trapped 
particles with and without mutual interactions [4] . First some preliminaries. 

Four lengths are involved in harmonic traps. They are: 

Li = (2nmuj 2 (3y 1/2 , L 2 = (h/muj) 1/2 , A = (27rh 2 (3/m)^ 2 , a, (1) 

where (3 = 1/kT, and L\ gives the order of the classical oscillation amplitude of a 
particle in the oscillator with energy kT, L 2 is the size of the ground state in the 
oscillator, A the thermal wave length and a the s-wave scattering length. For the 
recent experiments A >> \a\. We notice also that 

Li : L 2 : A = 1 : V2ne : 2ne, e = f3hw. (2) 

1. Thomas-Fermi method and the local density approximation (LDA), 
gaseous phase 

For small values of e and large number of particles, the Thomas-Fermi method 
for atoms can be adapted to the present problem: We divide space up into cells of 
volume larger than (L 2 ) 3 but smaller than (Li) 3 , and consider the potential energy 
V of the trap to be a constant in each cell. Each cell contains a collection of particles 
with mutual interaction described by a scattering length a. All particles in a cell have 
the same external potential V(r) per particle. If V — this problem for a cell has 
been studied [5,6] in the 1950s. For the present problem, in each cell we need only 
replace the fugacity z of the 1950 result by z exp(— f3V). 

This adaptation of the Thomas- Fermi method to the present problem has been 
used in a paper by Oliva (see Ref. [4]). We shall follow his terminology and call it the 
local density approximation (LDA). In this approximation, in the gaseous phase (i.e., 
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without BEC) in each cell, p(r) is given [5] by 

p(r) = \~ 3 g 3/2 (0[l-(4a/\)g 1/2 (()} (3) 

where 



oo 



C = zeM-PV), </„(*) = £rV (4) 
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and z is the fugacity. 

To order a/A we can also write (3) as 

p(r) = \- 3 g 3/2 (0 (5) 
i = z exp[-[3V(r) - 4a\ 2 p(r)}. (6) 

2. Validity of LDA for case a=0 

If a = 0, LDA [i.e., (5) and (6)] gives immediately an explicit expression for 
p(r) provided z < 1. If we pack in more particles, i.e., if we try to increase z beyond 
1, since £ cannot further increase at r = 0, BEC sets in at r = 0, in the same way as 
the original Einstein description of BEC in momentum space at p — 0. But this case 
of a = and V = ^muj 2 r 2 can be rigorously solved, allowing for an understanding 
of how LDA approaches the rigorous result. This is what we shall do in the present 
section. 

The density of atoms is p(r) = (r\D\r) where D is the density matrix equal 
to z e~P H {\ — z e _/3// ) _1 = J2T z e ~ ■ The matrix elements of e~^ Hl is explicitly 
known from [7]. (Using it we avoid the tedious process of summing over squares of 
Hermite polynomials.) We thus obtain 

oo 

p(r) = e 3/2 A- 3 ]T(sinh £e)- 3/2 z e exp[-a 2 tanh(fe/2)] (7) 
i 

where a = r/L 2 . The summand in (7) behaves at large £ like a geometrical series 
with the ratio of successive terms equal to 

Zl = z e ~ 3£/2 (8) 



Thus the summation converges at all r for z\ < 1 and becomes divergent at z\ — 1 
for all r. To study how it diverges we write the summation as Y^T a i where 

£n a t = £(£n z) - (3/2)£n(sinh is) - a 2 tanh(&r /2) 

= £[£n z - 3e/2] + (- £n2 - a 2 ) + C t 
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where Q -> as £ -> oo. Thus a* = 2 3 / 2 e - ^ + 2 3 / 2 e-^z^e^ - 1). 
Therefore the summation is equal to 

J2a e = 2 3 / 2 e- 2 ^^+ 2 3 / 2 e- CT2 5>f(e^-l). (9) 

1 — Z\ 

The first term on the right is the divergent part as z\ — > 1. We notice happily that 
it is exactly proportional to |-^o( r )| 2 where ipo is the normalized ground state wave 
function of the harmonic oscillator. Thus (7) becomes 

P(r) = -^\Mr)\ 2 + Pn(r) (10) 
1 - Zi 

showing that the BEC is in the ground state. p n is the normal fluid part of the density 
function p(r). It is the last term of (9) multiplied by A~ 3 £ 3 / 2 : 

p n (r) = A- 3 (2 £ ) 3 / 2 f>f j g— exp[- ( r 2 tanh( £ £/2)] - exp[-<r 2 ] j . (11) 

So far (10) and (11) are exact. Notice that for large £, the sum in (11) is convergent 
for zi — 1, unlike the first term of (10). Now we go to the case of £ << 1. The 
contribution to the sum in (11) for e£ > 1 is negligible. For e£ < 1, can drop the 
last term in the curly bracket and replace 1 — e~ 2fe by 2£e and tanh(e£/2) by e£/2, 
obtaining 

P(r) = -r^\Mr)\ 2 + \' Z 9^{z^ m ^). (12) 

1 — 21 

The second term here is the result of LDA with the replacement of z by z\. See (8). 
For e « 1, this replacement creates negligible errors. 

In the gaseous phase, i.e., 1 — zi~0(l), the second term in (12) dominates and 
the LDA is good. 
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In the BEC phase, the first term is the condensate, with 1 — zi^O(N~ 1 ). Thus 
we can put z\ — 1 in the second term, obtaining exactly the result of LDA, except 
for the fact that in the rigorous result the condensate has the spatial dependence 
of |-?/>o(V)| 2 — (const.) exp[— a 2 ] while in LDA the condensate has a delta-function 
dependence on f. But this is hardly surprising since in LDA each cell has a linear 
dimension large compared with L 2 , so that any structure of the order of L 2 is shrunk 
to a point. This fact also means that in LDA the cells should not be chosen to be < 
the order L 2 . 

To further check the error in (12) we evaluated numerically its error divided by 
the second term on the right for certain cases and found that the ratio is generally of 
the order of y/e. For example, if e — 0.05, the ratio is < 0.28, and for e = 0.01, the 
ratio is < 0.12. 

We summarize: As e — > 0, the LDA expression for p(r) approaches the rigorous 
result at every fixed a x = r/Li. 

We mention here that with a = 0, there is the well-known basic symmetry in 
the Hamiltonian between the coordinate and the momentum. To be precise, if we 
keep h, u and (3 unchanged, but replace m with {muj 2 ) 1 and switch x and p, the 
problem is unchanged (for a = 0). Thus we easily obtain the exact momentum space 
density distribution n(p). The condensate is of course in the state \ip (p)\ 2 . 

3. Condensate for the case a > in LDA 

For a > 0, the thermal equilibrium in each cell can be studied as in [6]. For 
large enough density, condensation takes place in the cell. The cell then consists of 
a saturated gaseous part with density po = \~ 3 g 3 /2(l) plus a super part with density 
p — p = p s . The free energy density then becomes, according to Eq.(33) of Ref. [6]: 

f(r) = -kT\- 3 g 3/2 {l) + 2a\ 2 kTp 2 - a\ 2 kTp 2 s . (13) 
From this free energy we obtain the chemical potential which should be equated 
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with kT£n[z exp(— (3V)\ giving 

kT In z-V = 47ia[p + p ]h 2 /m (r < r ). (14) 

At r = r — , p s = and p = p . Thus 

V(r) + 4nap(r)h 2 /m = V(r ) + 47vap h 2 /m. (14') 
This simple equation is valid in LDA for any trap potential V(r) and should be 

testable experimentally. 

Outside of r , (5) and (6) give the dependence of p on r. At r , £ = 1, and 
p = Po on both sides. I.e., p is continuous. The value of dp/dr is finite also but 
discontinuous at r , a fact already emphasized by 01iva[4]. So is d 2 p/dr 2 . 

4. Remarks 

The pseudopotential interaction [5] is , in q-number language: 

V = 4?ra ^ j- / dfidf 2 V' t (ri)V' t (r2)5(fi - f 2 )-^-[ri 2 ^(ri)^(f 2 )]. (15) 
m 2 J ar 1 2 

The commonly used delta-function interaction [5] is, also in q-number language, 

Anah 2 1 



J df^\r^^(r)^(r)^(f). (16) 



m 2 

(a) We point out that there has been great confusion in the theoretical literature of 
factors of two about V pp and Vs. Furthermore the approximate one-particle c-number 
equation written down from (16) also has confusions of factors of two in the literature. 

(b) Strictly speaking (16) does not make sense: The two-body interaction A8{f\ — f 2 ) 
for A < is not defined and for A > is equivalent to zero. To prove the first 
statement, we calculate the s-wave phase-shift for such a case. The interaction is 
considered as the limit R — > of a potential well of magnitude U and radius R, 
keeping UR 3 = negative constant. The wave function in the center-of-mass system at 
R— is R' 1 sin[(— mU/h 2 ) 1 ^ 2 R]. The argument of sine approaches infinity as R — > 0. 
To prove the second statement we do the same calculation and find all phase- shifts 
to be zero. 
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(c) Although (16) is strictly meaningless, the perturbation calculations based on (15) 
and (16) in Refs.[5, 6], give the same result to the order a 1 , and is meaningful to that 
order. The considerations of the present paper are based on these results to order a 
and are thus meaningful. 

How about higher orders? This is a matter of some subtlety. It was found in 
Ref. [5] that to the second order in a, (16) gives infinity, confirming its sickness, but 
that (15) gives a meaningful answer which is proportional to a/L where L is the size 
of the box. To order a 3 (15) gives an energy containing a term equal to (constant) 
times N 3 a 3 /L 5 which approaches infinity as L — > oo while p = N/L 3 = constant. 
This divergence was later removed [8] by a method of summation over most divergent 
terms in a perturbation expansion: A summation starting with this divergent term 
leads to the convergent expression (AixaNp) (128/15) (^p/n) 1 ^ 2 . This result was later 
confirmed by several authors who extended it to even higher orders. 

The method of Ref. [8] was further extended to cover finite temperatures [9]. 
These and related developments had been summarized in a brief talk [10]. 

(d) It is tempting to interpret (5) and (6) as indicating an additional effective potential 
V a = kTAa\ 2 p = 8iraph 2 /m, and to interpret (14') as indicating an additional poten- 
tial of V a = Anaph 2 /m. Such interpretations must be used with care. In particular, 
it is not correct to assert that the condensate is in the ground state of V + V a . 

The work of CNY is supported in part by an NSF Grant PHY-9309888. The 
work of LHY is performed under the auspices of US DOE. 
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FIGURE CAPTION 

Fig. 1 Example of p(r) as a function of r for a > according to Eq. (14'). The 

curve is for a harmonic trap. Notice that the first and second derivatives are both 
finite but discontinuous at r . po = A~ 3 g 3 / 2 (1) is the normal part of p, and p s the 
condensate part. For smaller total number of particles (= A), r becomes smaller. It 
eventually shrinks to zero and the condensate disappears. 
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